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Abstract 

In graphene, after tlie electric field is turned-on, the ballistic acceleration of charge carriers 
up to the monochromatic optic phonon energy generates a back-and-forth motion of the whole 
distribution function between the zero point energy and the phonon energy. This effect is 
predicted to manifest in damped terahertz oscillations of the carrier drift velocity and average 
energy. The most dramatic feature of this transient phenomenon is the onset of momentum-free 
areas suiTounded by high momentum probability in phase space, and smooth steps or teiTaces 
in the distribution function. This dynamical effect that only takes place within a voltage and 
sample length window, is the direct consequence of the interplay between the electric force 
and the randomizing nature of deformation potential optic phonons in the linear band structure 
of graphene. 
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Since its isolation in monolayer sheets, graphene has emerged as a new and ideal two-dimensional 
(2D) material with exotic electronic properties.^ The particular nature of its band structure with a 
linear dispersion between energy and momentum E = hvf\k\, where Vf ~ lO^cm/s is the Fermi 
velocity of electron and holes, results in a zero-gap material, where the overlap between the va- 
lence and conduction band is reduced to a single point around the K and K' symmetry points of the 
Brillouin zone.^' " Hence, charge carrier dynamics can be described by a formalism similar to Dirac 
relativistic equation, where the speed of light is replaced by Vf which provides new opportunity 
to observe "quantum electro-dynamic" effects in nanoscale systems. For electronics applications, 
graphene is similarly attractive because of the linear energy-momentum dispersion, for which all 
charge carriers move with the velocity, Vf much larger than in conventional semiconductor mate- 
rials,^ which anticipates faster time response. 

Indeed, the carrier mean free path for collisions with low energy acoustic phonons (AP) can 
reach several micrometers, while efficient scattering with monochromatic optic phonons (OP) 
occurs at much larger energy {hcOop ~ 0.2eV) than in conventional semiconductors {hcOop ~ 
0.04ey).^^ Therefore, in electric fields F high enough (to escape AP or low energy scattering), 
charge carriers can experience quasi-ballistic run-away until they scattered with efficient OP emis- 
sion once E > hcOop.^ The coherent aspect of this process i.e. the coherent acceleration of the 
carrier distribution function followed by quasi-instantenous carrier relaxation by OPs at high en- 
ergy is expected to produce oscillations of the carrier velocity with a periodicity given by T^^ = 
hcOop/eFvF.^ 

In this paper we investigate carrier dynamics in graphene single layers, once the electric field 
has been tumed-on, to determine the conditions of occurrence of current oscillations. Indeed, 
the thermal broadening of the initial distribution introduces a decoherence in the OP relaxation 
amongst carriers, which asides from low energy scattering, causes inherent damping of the os- 
cillations, especially at high fields, where carriers overshoot the OP energy before emission. For 
this purpose, we use the Boltzmann formalism, and solve for the time-varying carrier distribution 
in the presence of OP scattering. ^"^'^ We also account for low energy scattering such as impu- 
rities and acoustic phonons, by using the relaxation time approximation. In weak concentrations 
{uc < lO^^c/M^^), electron-electron interactions do not play a major role in transport in graphene, 
and are not included in this analysis. We provide an analysis of the interplay between applied elec- 
tric fields and strength of the low energy scattering rates for the onset of oscillations. 

Model. Let us consider a system of electrons in the graphene conduction band under the influ- 
ence of an electric field along the x-direction. We divide the momentum space into two regions, I 
{k < kc) and II {k > kc) separated by a circle of critical momentum kc = (Oop/vp, that corresponds 
to the electron kinetic energy equal to hcOop (Fig. la). 

In region I electrons undergo quasi-ballistic acceleration and weak scattering by low energy 
mechanisms (e.g. impurities or APs) until they reach region II where they lose their energy by 
OP emission, and scatter back to region I. In our model, the electric field is assumed to be low 
enough so that electrons are scattered efficiently from region II to region I by OP emission with 
little probability to reach E > IhcOop. 

In the two regions, the time-dependent Boltzmann equations read: 
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Figure 1: a) Schematics of caiTier quasi -ballistic acceleration and OP scattering in 2D k space. The circle 
of radius kc is the locus of all points in k space con^esponding to the earner energy fiCOop. b) Schematics of 
electrons scattered by OP's from two different positions of the distribution function in region II. The dotted 
circles in region I represent the two areas of high momentum probability where electrons are more likely to 
land. 
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k' 

where fi(kj) and fji(k,t) are the time-dependent momentum distribution functions in regions I 
and II respectively, and F is the electric field applied.^ The first terms in RHS of equations (1) 
account for low energy scattering mechanisms where, 

Mk) = (2) 
l+exp(^) 

is the Fermi-Dirac equilibrium distribution function, and T is the relaxation time (for the sake of 
simplicity, we assume T is ^-independent and we vary it's value compared to the OP scattering 
rate). S{k' Jc) is the OP transition rate from a state with momentum k' to the state with momentum 
~k, given by " - 

..r r/^ TtDo'{Nq+\)d{E'-E + h(Dop) 

S{k,k) = (3) 

oAoop 

where Do is the optical deformation potential, a is the mass of the graphene sheet per unit area, 
and A is the area of the sheet. 

The second term on the right-hand side (RHS) of Eq. lb is the carrier depopulation by OP 
emission, while in Eq. la it is the corresponding carrier repopulation at low energy. Here we 
neglect OP absorption processes since hcOop ^ kgT , even at room temperature where the phonon 
occupation number A^^ is negligible, so OPs only scatter electrons from region II to region I. 

We set the initial distribution {t = 0) as the Fermi-Dirac distribution, 

f,{k,,ky,t=0)=Mk,X) (4) 

and solve Eqs.l iteratively by noticing that for times <t < hkc/eF, the inside distribution drifts 
towards the critical circle with "speed" ef/h. During this trip, one assumes that all carriers are 
in region I, corresponding to /// ^ 0. Eq. la is just a differential equation for which a solution is 
readily obtained. Then we solve Eq. lb for k > kc, which is of the same form as Eq. la but with the 
additional OP emission term. Here we use the boundary condition fi{k = kc,t) = fii{k = kcj). ' ' 
For t > hkc/ eF, we substitute the solution fnik, t) into the integral of the right hand side of Eq. la 
to start the same procedure for later times. 

As the electron population moves back and forth between regions I and II, each time undergoing 
more dephasing and broadening due to the finite duration of the OP emission process, we write 
the distribution function in each region at any time t corresponding to the n''^ trip toward the 

critical circle, i.e. for n = integer[eFt j tike] -l- 1 , as a superposition of distributions ff' {k^, ky, t) and 

fli\k^,ky,t) of individual/^'' "trip" , 

Mk,,kyj) = f^fi'\k,,kyj) (5a) 
(=1 
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and 

fn{k,,ky,t) = f^fl'l^{k,,ky,t) (5b) 

i=l 

We note that the concentration ric in the conduction band is found to be ~ 8 x lO^^cm"^ for 
Ef = (at the Dirac point) used throughout this analysis. For Fermi levels at Ef = kgT and 
Ef = IkgT above the Dirac point, the carrier concentrations would be tic ~ 1.8 x lO^^cm^^, and 
He ~ 3.4 X I0^^cm~^, respectively. 



a) b) 




Figure 2: a) 3D plots of the distribution function at different times ([5 = 0, 1, 1.5, 2, and 2.5) for F = 
lOOOV /cm and vanishing low energy scattering (7 = 0). b) CoiTcsponding depopulation and repopulation 
rates for the same times as in (a) 

Results. We introduce the dimensionless time parameter /3 given by 

^ ~'t* 

where 

ef 

is the approximate time taken for the distribution to complete one "trip". 

Fig. 2a shows the time evolution of the distribution function distribution (DF) for F = lOOOV / cm 
(t* ^ 2ps), in the absence of low energy scattering (7 = To/ T = 0). At /3 = 0(r = 0) the initial dis- 
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tribution is the tail of the Fermi distribution in the conduction band (Eq. 12), which is Maxwellian- 
like, and drifts along the direction (opposite to the F-field) to reach the OP energy at /3 = 1 
(second panel). At that time, a "hump" appears around k = (not yet visible on the DF graph) as 
the front electrons reach region II, where they emit OP's and scatter back inside region I, as seen 
in the second panel of fig. 2b. As time progresses /3 = 1.5, the "hump" develops into a dimple, 
while the remaining high energy electrons from the first trip continue their drift in region II, where 
they experience strong OP depopulation. As seen in the third panel of fig. 2b, the corresponding 
repopulation rate at low energy exhibits a crater-like shape, which is due to the randomizing na- 
ture of the deformation potential OP scattering. Indeed, as schematically shown in fig. lb, after OP 
emission by high energy electrons, all \k — kc \ -values become equiprobable, which forms a drifting 
circle of high repopulation rate, with increasing k-radius as first trip electrons penetrate deeper in 
region II. This crater-like feature of the repopulation rate is the primary cause of the dimple in the 
low energy distribution function, which are areas in k-space where electrons have low probability 
to scatter. As time progresses, the DF "dimple" evolves into a crater-like shape, (fig. 2. a, fourth and 
fifth panel), and the successive depopulation-repopulation OP processes overlapping at low energy 
with different amplitudes form also smooth terraces in the low energy tail of the distribution as it 
approaches steady state (fig. 2. a , fourth panel). These morphological effects in electron distribution 
are unique to graphene as a result its linear band structure and the interplay of the quasi-ballistic 
acceleration and relaxation by high energy monochromatic OP's. 




Figure 3: 3D plots of the distribution functions at j8 = 2.5 for 7 = 0, 0.05, 0.1, and 0.5. The applied field is 
F = \OmV/cm. 
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Fig. 3 shows snapshots of the distribution function at j8 = 2.5 for varying low energy relaxation 
time, expressed in terms of 7 = Tq/t. From the figure, one can see that the crater in the DF that 
occurs for 7 = progressively disappears as 7 increases. Indeed low energy scattering in region I 
re-distributes charge carrier momenta around k = 0, especially in the crater center. For this reasons, 
the amplitude of the distribution in region I increases around ^ = 0. Also, the distribution amplitude 
decreases in region II (k > k^) as carriers spend more time in region \ {k < kc). Similarly the 
distribution recovers a streaming profile along the ^^-direction as conventional semiconductors. ' 
Nevertheless, even for strong damping (7 = 0.5), the distribution is characterized by a jagged 
profile, which contains the front and back ridges of the crater remnant still caused by the cumulative 
effects of the OP scattering for backward and forward carrier relaxation. 




Figure 4: a) Current density as a function of time for different values of the low energy scattering parameters 
7. Solid (7 = 0), Circles (7 = 0.05), Crosses (7 = 0.1), and Diamonds (7 = 0.5). The applied field is 
F = lOOOV /cm. Inset: Cross-sections of the distribution at jS = 3 and kx = 0.5kc for the corresponding 
values of 7. Arrows; initial distribution function at k^ = 0. b) Current density as a function of time for 
different fields and 7 = 0.1: Solid (F = 1500V /cm), Circles (F = lOOOV /cm), Crosses (F = 500V /cm). 
Inset: Cross-section of the distribution function at /3 = 3 and kx = 0.5kc for the corresponding values of F. 
Diamonds; initial distribution function at kx = 0. 

The current density on the plane is given by, 

7,(0 = -evF £/ A t) cos((^)) (6) 
k 
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where ^ is the angle between k and the axis. Fig. 4a shows the current density as a function 
of time (/3) for F = lOOOV/cm and different low energy scattering rates. For weak scattering, the 
current density overshoots its steady state value through damped oscillations, as a result of the 
back and forth motion of the distribution function (fig.l). The weaker the scattering, the higher the 
current overshoot. For strong low energy scattering, the current converges monotonically toward 
its steady state value without oscillations. ^ Notice the stronger the scattering, the lower the steady 
state current value. The insert shows the corresponding DF cross-sections at kj^ = 0.5 and /3 = 3 
relative to the initial distribution at fc^ = 0. As low energy scattering increases, the DF becomes 
narrower, taking a "streaming" profile in the electric field direction. This is due to the fact that 
as the low energy scattering increases, fewer electrons reach the high energy ( E > ficOop) region 
II, reducing the number of electrons scattered back to the low energy region I, which narrows the 
distribution. 

Fig. 4b shows the current density as a function of time for three different field values, with 
7 = 0.05 in each case. One observes that the oscillation period Xgr = ficOop/eFvf scales with the 
inverse of the electric field F . Also quite expectedly, the overshoot value increases with electric 
fields, but the damping is also enhanced with electric fields, which is due to the fact that the 
electron distribution penetrates the high energy(£' > hcOop) region(II) farther than hcOop, stretching 
the carrier relaxation by OP emission, which in turn broadens the DF in region I along the field, 
thereby reaching steady state quicker. The insert shows the DF cross-sections at k^ = 0.5 and /3 = 3 
for the different fields, relative to the initial distribution at k^ = 0. As the field is increased, the 
distribution function broadens, because more electrons reach the high-energy {E > fiOJop ) region 
II, and scatter back into region I broadening the distribution in the process. 

The average energy density of carriers in the conduction band is given by, 

E{t)=hvFY,kfCk.t) (7) 
k 

where f{k,t) is the overall normalized distribution function, and the summation is taken over all k 
in regions I and II. At /3 = 0, the average energy per charge carrier is found to be E(0) ^ l.lkgT . 
This value is roughly twice larger than ksT expected for two-dimensional systems, and is the direct 
consequence of the linear energy dispersion in graphene, by contrast to the parabolic dispersion 
in normal 2D systems. Fig. 5a displays the ratio £ = £'(/3) /E{0) for F = \kV /cm and different 
low energy scattering rate (7). As expected, the energy converges to higher values as low energy 
scattering is reduced. In addition, the convergence toward steady state occurs through damped 
oscillations, even for 7 = as a consequence of the back and forth motion of the DF between the 
OP energy and the carrier zero-point energy. Quite clearly the oscillation period is given by t* for 
all 7. It is also seen that the oscillations persist even with significant low energy scattering. Fig. 5b 
displays the normalized energy £ as a function of time for different fields and 7= 0. 1 . As expected, 
carrier energies reach higher values as the field is increased, and the oscillation period decreases 
(larger j8 -period). 

Discussions. We have provided a transient analysis of the onset of current oscillations at the 
electric field turn-on caused by the back and forth motion of carrier distribution function between 
the zero-point energy and OPs in the presence of varying damping mechanisms. In this context we 
point out the anomalous shape of the carrier distribution as an interplay between ballistic accel- 
eration and deformation potential OP emission in the transient regime. If OP-limited current os- 
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Figure 5: a) Normalized value of the carrier energy as a function of time for different values of the low 
energy scattering pai^ameters 7. Solid (7 = 0), Circles (7 = 0.05), Crosses (7 = 0.1), and Diamonds (7 = 
0.5). The applied field is F = lOOOV /an. b) Same but for different applied fields and 7 = 0.1: Solid 
{F = 1500V /cm), Circles (f = lOOOV /cm), Crosses {F = 500V /cm). 
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cillations have been predicted in GaAs, and indirectly observed in slightly n-doped InSb, their 
manifestation in graphene is different in several respects: First, owing to the linear carrier energy- 
momentum dispersion, the oscillation periodicity is given by Tg^ = hcofp/eFvp in graphene, while 



in GaAs parabolic conduction band it is expressed as XgoAs = 2m*hco2p^^ / , which neverthe- 
less yields similar values, since the small effective mass, and the OP phonon frequency in GaAs 
compensate for the large v/r- ' Second, in III-V semiconductors, the OP polar nature focuses 
the low-energy repopulation along the axis, which provides a "streaming" profile to the carrier 
distribution instead of a crater-like shape in this case. Finally, in compound semiconductors, the 
oscillation onset is restricted by two conflicting conditions: On the one hand, the low value of the 
OP energy {fiCDop '■ O.OAeV) is comparable to the thermal broadening of the carrier distribution at 
room temperature so that the back and forth motion of the distribution between the optic phonon 
and the zero point energy is immediately damped. On the other hand, at low temperature, 
ionized impurity scattering becomes dominant and produces strong damping which can only be 
reduced by lowering the dopant density, thereby lowering the carrier density, and weakening the 
oscillation amplitude. In this respect, the high conductance of graphene, and the high optic phonon 
energy provide the conditions for room temperature observation. 

Low energy scattering should however still be minimized. Usually, for experimental studies 
and device applications, graphene-layers rest on a dielectric substrate, or are confined between 
two dielectric slabs. In this case the presence of high K-dielectrics sandwiching the strictly 
2-D graphene-layer may be used to screen charged impurity that may reduce low energy elastic 
scattering. ' However, dielectrics also contain interface and remote static charges that may offset 
dielectric screening."^ Moreover, the interaction between the 2D carriers in graphene and low 
energy Remote Interface Phonon (RIP) arising from the proximity of the substrate introduces 
new scattering sources. ^^'^^ Therefore in general, suspended graphene-layers avoiding the RIP 
influence may be preferable. - 

In this context to be observable at room temperature, the velocity oscillations should also take 
place within a parameter window. On the one hand the process requires eV > fiCOop where V is the 
external bias, to reach the OP energy, and on the other hand, the field should be large enough for 
carriers to escape low energy scattering, but hcOop/ eFvfXop ^ 1 where 1 / top is the OP scattering 
rate, so that charge carriers don't penetrate the high-energy region E > ficOop, scatter immediately 
after they reach the OP energy, which maintains the coherence of the distribution function. These 
requirements impose a lower and upper bound on the electric field, i.e. O.SkV /cm < F <^ 5kV /cm 
(in graphene), and a lower bound on the sample length L > VfTop (> l/^m for Top < Ips), but L 
should be smaller than a few values of A = hcOop/eF, so as to prevent oscillation damping. 

One important issue for the validity of this analysis is the effect of leakage current due to 
electrons in the valence band crossing to the conduction of carriers through the Dirac point. Indeed 
it has been shown that there is still a minimum conductance (G ~ 4e^/h) between the two bands 
despite the singular nature of the Dirac point. This value was later measured to be G ~ e-^/h. ' 
However, this effect becomes important only when the graphene layer width w > 23/im, for which 
the band gap Eg < OASmeV are vanishing. Moreover, thermal effects such as acoustic phonon 
absorption by carriers in the valence band, are forbidden by conservation of both, energy and 
momentum. As for OP absorption we have shown earlier their occupation number is also negligible 
over the time scale considered in this analysis. 
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Supporting Information 
Solution of Eqs. la &; b. With the initial distribution [t = 0) given 
the Fermi-Dirac distribution, 

fi{kx, ky, t = 0) = fo{kx, ky) 

the solution for the "first trip" (0 < t < hkc/eF) in region I [k < kc) is, 

fj^\kx,ky,t) = —exp{-'^^^) f dzexp{—)fo{z,ky) 



T 



+ fo{kx--,ky)exp{--) (S.2) 

Uq t 

with 

h 

= 

eF 

where I/z/q is the "speed" in /c-space with which the distribution drifts in the 
field towards the critical circle. 

High Energy Distribution Function {k > k^). Since hvfk ^ ksT 

for k > kc even at room temperature, the approximation fo{k) ~ will be 

used for E > hwop- By setting a new variable u = t — hk^/eF, Boltzmann 

equation (l.b) in region II reads, 

dgii{k-u) - Jl 1 ) . 

- Uogiiik; u}<-^ — - > (S.3) 



dkx [t Top{k) 

where 

gn{k\ u) = gn{kx, ky\ u) = fn{kx, ky] u + Vokx) 
and the total OP emmission rate in Eq. lb {k > kc) is. 



Top{k) ^ ' 2ahwopVf To [kc 

1 



If we use hbUop = 0.2eV, = l^eV/°A, a = 7.61 x IQ-'^kg/m'^, then ^ 
0.32psec. The solution of Eq. S.3 is then, 

gn{h: ky-, u) = gii{k°, ky\ u) exp(-r(A;^, ky)) (S.5) 

where 



T{k^, ky) = 1^0 [ dk^\--\ \-- \ 



and 



kx — \/kc ky 



The solution (5". 5) is valid only in the region —kc < ky < kc, which is the 
region of interest. Going back to the original variables, the general solution 
reads, 

fii{kx: ky, t) = fn{kx: ky, t - u„{kx - k^) exp(-r(A;^, ky)) (S.6) 

The above solution describes a distribution function that drifts with "speed" 
eF/h'm. momentum space in the direction of the electric field, while decaying 
exponentially by OP emission. 

Low Energy Distribution Function {k < kc). By using the same 
substitution as in the previous paragraph, the Boltzmann equation (l.a) 
inside the critical circle (region I) reads, 

dgi{k] u) i^o{9i{k] u) - fo{k)} 



dk. 



+ Z{k;u) (S.7) 



X 



where 



Z{k; u) = ^° + f d(j)'fii{k + kc, 0', u + Vokx 

A-KGUeFvp J-a(k) 



2 



and 



k 

a{k) — arcsm{-- — ^-—-} 



The solution to (5". 7) is, 

gi{k;u) ^ ex-p{-'^^) I d2;exp(— ) i ^°f"^^'^y> + Z{z,ky;u)\ 
^ J-k° r I r ) 

+ n{ky;u)exp{-'^) (S.8) 

T 

in which the function fl{ky;u) is determined by initial conditions. Going 
back to the initial variables, 

//(^x,^?/,*) = exp(-^^^) f dz exjp{—) + Z{z,ky,t - Uokxrl 

J-k° T y T ) 

+ ^{ky,t- Uokx)exp{-'^^) (S.9) 

r 

By using the matching condition, 

fi{klM = fn{kl,ky,t) (S.IO) 

and initial conditions, the solution for the distribution function in regions I 
and II is readily obtained. 

Total Distribution Function. Finally, we write the distribution func- 
tion in each region at any time t corresponding to the n*'* trip toward the 

critical circle, i.e. for n = integer[eFt / hkc] + 1, as a superposition of distri- 
butions ff\kx,ky,t) and fj^j{kx,ky,t) of individual i*^ "trip" , 

n 

fi{kx,ky,t) = ^ ] /j \kx, ky, t) (S.lla) 

i=l 
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and 

n 

fii{h,ky,t) = Y,fihk.,ky,t) (S.llb) 
i=l 

The normalization condition 

= x{t) I Yl /^(^'' + E /^^(^'' ^) f (S.12) 

determines as Eqs. (S*.!!) do not exactly conserve the number of carriers 
fic- Indeed there is a decrease in the number of carriers in the system because 
at high fields electrons run away beyond region II i.e > 2f}UJop without being 
scattered. Another reason is that the summation in Eqs. (-S'.ll) is hmited 
to the n = 3 terms thereby ignoring the contribution of later "trips" (n > 4 
terms) which might be relevant at higher fields considered in the analysis. 
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